/*=================================================================================
 * Class Reaction:	Class for reaction process
 *					Read file -> Construction -> Initialization -> Reaction
 *=================================================================================*/

#ifndef REACTION_H
#define REACTION_H

/*---------------------------------------------------------------------------------
 * Standard C++ library headers
 *---------------------------------------------------------------------------------*/

#include <fstream>
#include <iostream>
#include <string>

/*---------------------------------------------------------------------------------
 * Local headers
 *---------------------------------------------------------------------------------*/

#include "Array.h"

class Reaction {
private:
	/* x axis */
	int ni;								/* Number of nodes */
	Array<double> xnode;				/* The mesh nodes */

	/* y axis */
	int nj;								/* Number of nodes */
	Array<double> ynode;				/* The mesh nodes */

	int bc;								/* The number of each boundary grid */

	int NS;								/* Number of species */
	int NR;								/* Number of reactions */
	int TB;								/* Number of three-body reactions */

	Array<double> Stoi_F;				/* Stoichiometric coefficient of reactants */
	Array<double> Stoi_B;				/* Stoichiometric coefficient of products */
	
	/* Forward reactions */
	Array<double> Af;					/* Pre-exponent factor */
	Array<double> Bf;					/* Temperature index */
	Array<double> Eaf;					/* Activation energy */
	
	Array<double> React_TB;				/* Coefficient of three-body reactions */

	Array<double> Coeff0;				/* NASA Thermochemical polynomial coefficient (T<1000) */
	Array<double> Coeff1;				/* NASA Thermochemical polynomial coefficient (T>1000) */
	
	const double R = 8.31434;			/* gas constant of ideal gas J/(mol*K) */
	const double Ru = 1.987;			/* gas constant of ideal gas cal/(mol*K) */
	const double P0 = 101325;			/* one standard atmosphere */

	Array<double> Mw;					/* Component molecular weight */
	Array<double> Ri;					/* Component gas constant */

	Array<double> Mc;					/* Component mole concentration */
	Array<double> Mr;					/* Component mole ratio */
	Array<double> Mi;					/* Component mole fraction */
	Array<double> Yi;					/* Component mass fraction */
	Array<double> Di;					/* Component density fraction */

	double T_local;						/* Local temperature */

	Array<double> Hi;					/* Enthalpy of each component */
	Array<double> Si;					/* Entropy of each component */
	Array<double> Gi;					/* Gibbs free energy of each component */

	Array<double> KF;					/* Forward reaction rate constant */
	Array<double> KB;					/* Reverse reaction rate constant */
	Array<double> Kp;					/* Chemical reaction equilibrium constant expressed in terms of pressure */
	Array<double> Kc;					/* Chemical reaction equilibrium constant expressed in terms of concentration */

	Array<double> RR_F;					/* Forward reaction rate */
	Array<double> RR_B;					/* Reverse reaction rate */
	Array<double> R_TB;					/* Correction coefficient of three-body reactions */
	Array<double> RR;					/* Net reaction rate */

	Array<double> Wi;					/* Component mass production rate */
	Array<double> CMS;					/* Component mass production rate */
	Array<double> WJH1, WJH2;			/* Intermediate variables */
	Array<double> MD;					/* Diagonal of the Jacobi matrix */

	Array<double> P, Q;					/* The variables for Trapezoid formula */

	Array<double> Nchem;				/* The number of chemical iteration step of every mesh node */
	int NchemNow = 1;					/* The current maximum Nchem */

public:
	friend class Euler;					/* Friend class -> access to use private object */

	/* Constructor */
	Reaction() = default;

	/* Read the reaction model file and thermo properties */
	void ReactionRead(char* reaction_model, char* thermofile, Array<double> xnode, Array<double> ynode, int bc);

	/* Construction */
	void ReactionConstruction();

	/* Calculate some initial parameters */
	void ReactionInitial();

	/* Reaction process using Trapezoid formula */
	Array<double> Trapezoid(Array<double> Mc, Array<double> Di, Array<double> T, double dt);
	Array<double> Trapezoid(Array<double> Mc, Array<double> Di, Array<double> Yi, Array<double> T, double dt, int step);	/* Adaptive */

	/* Reaction process using IMEX with diagonalized matrix */
	void Diagonalized(Array<double> Mc, Array<double> Di, Array<double> T, Array<double> Partial_T);

	/* Reaction process using IMEX with the exact matrix */
	void FullMatrix(Array<double> Mc, Array<double> Di, Array<double> T, Array<double> Partial_T);
	
	/* Reaction process using Cantera or DNN */
	Array<double> ReactionS(int type, Array<double> Di, Array<double> T, Array<double> P, double dt);

	/* Destructor */
	~Reaction() { ; };

	/* Get the specific heat at constant pressure */
	double GetCpi(double T, double R, int SP, Array<double> Coeff0, Array<double> Coeff1);

	/* Get the enthalpy */
	double GetHi(double T, double R, int SP, Array<double> Coeff0, Array<double> Coeff1);

	/* Get the entropy */
	double GetSi(double T, double R, int SP, Array<double> Coeff0, Array<double> Coeff1);

	/* Get the component mass production rate */
	Array<double> GetWi();

	/* Get the diagonal of the Jacobi matrix */
	Array<double> GetMD();

};
#endif
